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Abstract 

We propose an efficient method to evaluate callable and putable bonds under a wide class of interest 
rate models, including the popular short rate diffusion models, as well as their time changed versions with 
jumps. The method is based on the eigenfunction expansion of the pricing operator. Given the set of call 
and put dates, the callable and putable bond pricing function is the value function of a stochastic game 
with stopping times. Under some technical conditions, it is shown to have an eigenfunction expansion 
in eigenfunctions of the pricing operator with the expansion coefficients determined through a backward 
recursion. For popular short rate diffusion models, such as CIR, Vasicek, 3/2, the method is orders of 
magnitude faster than the alternative approaches in the literature. In contrast to the alternative ap- 
proaches in the literature that have so far been limited to diffusions, the method is equally applicable to 
short rate jump-diffusion and pure jump models constructed from diffusion models by Bochner's subor- 
dination with a Levy subordinator. 
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1 Introduction 

A large fraction of all corporate and sovereign bond issues in the global financial markets have embedded 
options. The call option allows the bond issuer, such as a corporation or a government, to buy the bond back 
from the bond holder (call the bond) for pre-specified call prices at some pre-specified times prior to bond's 
maturity. This allows the bond issuer to refinance the bond if interest rates decline. The put option allows the 
bond holder to sell (put) the bond back to the bond issuer for pre-specified put prices at some pre-specified 
times prior to maturity. This allows the bond holder to re-invest the proceeds into a bond with a higher 
coupon if interest rates rise. The bond with both a call and a put option can be analyzed as an insta nce of a 
stocha stic game with stopping times (also known as Dynkin games as they have been introduce bv iDvnkin 



(|1969f ) as a generalization of optimal stopping problems) driven by the underlying stochastic interest rate 
model. The bond issuer and the bond holder are opposing players whose opposing optimal strategies are to 
minimize and to maximize the bond value, respectively. When the call and put decisions can be made at 
discrete times (typically an advance notice has to be given to the other party ahead of a coupon payment 
date, when the option is exercised), this sets up a stochastic game with stopping times in discrete time. The 
value function and the optimal call and put policies can then be determined by solving Bellman's dynamic 
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programming backward induction, starting from maturity and rolling back recursively through the decision 
times. Developing solution methods for this problem is of significant practical importance. 

The pricing of bonds wi th embedded options h a s attr acted considerable interest in the literature over 
the years. It goes back to iBrennan and Schwartz! (1977) who used a finite-difference approach in time- 
homogeneous diffusion models. Later iButtlerl ( 19951) showed that when finite-difference methods were used 
to price callable bonds under the Vasicek model, the presence of slowly deca ying oscillations in the s o lution 
after each coupon/call date resulted in poor numerical accuracy. This led iButtler and Waldvoeell (jl996l ) 
(BW) to develop an alternative method for pricing callable bonds under the Vasicek and CIR models utiliz- 
ing the explicit form of the Green's function in these models. The method relies on the interpolation of the 
value function and on the numerical qu adrature procedure for the integration involving the value function 
and the Green's function. More recently Id'Halluin et al. ( 200ll ) (DFVL) showed that finite-difference meth- 
ods for these problems can be stabilized via van Leer flux limiter, appropriate non-uniform time stepping 
schemes, and careful consideration of boundary conditions and presented numerical experiments demon- 
strating that properly formulated finite-difference methods, in fact, outperformed alternative approaches. 
Other recent works on applying numerical PDE meth ods to callable a nd p utable bonds in diffusion intcr- 
est rate models include Farto and Vazque"3 ( 20051 ) and de Frutosl (2008). In Farto and Vazque"3 (2005), the 
convection dominated diffusion equation is solved numerical ly by combining the characteristics method with 
piecewise-linear Lagrange finite elements, de Frutosl (|2008l ) (F) proposes a spectral numerical method for 
pricing callable bonds. The holding value function is approximated as a finite summation involving the La- 
guerre polynomials, and the problem is converted to a stiff sys tem of ordinary di f ferent ial equations (ODEs) 
for the time-dependent coefficients of the Laguerre expansion. Ben-Ameur et al. ( 20071 ) (BBKL) propose an 
alternative dynamic programming approach not based on PDEs in which a piecewise linear approximation 
is used for the value function and the exact transition probability is used to compute the expectation of the 
discounted piecewise linear approximation of the value function in the cases of CIR and Vasicek short rate 
models where the exact analytical expressions are available. All of these papers provide numerical experi- 
ments illustrating computational performance of their respective methods on the same example of a Swiss 
callable bond. This provides a natural comparison benchmark. 

The present paper proposes an efficient method to evaluate callable and putable bonds under a wider- 
class of interest rate models than any of the previous approaches, including the popular short rate diffusion 
models, as well as their time changed versions with jumps, including both jump-diffusion and pure jump 
models with state-dependent jumps. The method is based on the eigenfunction expansion of the pricing 
operator. We show that, under some technical conditions, the callable and putable bond pricing function 
has an eigenfunction expansion in eigenfunctions of the pricing operator with the expansion coefficients 
determined t hrough a backwa rd recursi on. For popula r shor t rate diffusion mode ls, such as Cox-Ingersoll- 



Ross (CIR) (|Cox et all 119851 ) . Vasicek (jVasice: 



: p opula r 



3/2 dAhn and Gaol . Il999l) . we demonstrate on the 



test case of the Swiss bond used in the previous studies that the method is orders of magnitude faster than 
the alternative approaches in the literature. In addition, in contrast to the alternative approaches in the 
literature that have so far been limited to diffusions, the method is equally applicable to short rate jump- 
diffusion and pure jump models constructed from diffusion models by Bochner's subordination with a Levy 
subordinator. 

The strength of the eigenfunction expansion method is that the value function is constructed globally in 
state and time with no need for discretization of either state or time variables. The only approximations 
involved in the computation scheme based on the method are the truncation of the infinite eigenfunction 
expansion (that, under some technical conditions, is uniformly convergent with uniformly controlled trun- 
cation error) and the numerical solution of a non-linear equation to determine the stopping boundary at 
each step of the backward recursion solved by the fast-converging bisection algorithm. Another strength of 
the method is that it can be seamlessly applied to both jump-diffusion and pure jump interest rate mod- 
els obtained from diffusion models by subordination. Semi-analytical methods of BW and BBKL are not 
suitable to handling jump-diffusion and pure jump models since no analytical solutions are available for 
transition probabilities and Green's functions in these models. Numerical PDE methods, such as finite- 
difference and finite element methods, can, in principle, be extended to handle partial integro-differential 
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equations (PIDE) arising in jump-diffusion and pure jump models, but at substantial costs both in the 
theoretical complexity and structure of the schemes and in their computational implementation and com- 
putational performance. Moreover, most of the existing instances of applying numerical PIDE methods in 
computational finance have been limited to Levy processes with state-independent jumps. In contrast, the 
eigenfunction expansion method is capable of handling models with state-dependent jumps, such as mean- 
reverting jumps in the interest rate. Both at the theoretical and computational level, moving from a pure 
diffusion model to a jump-diffusion or a pure jump model obtained from the diffusion model by the time 
change with respect to a Levy subordinator amounts to no more than replacing the eigenvalues e _A ™* of 
the pricing operator in the original diffusion model with the eigenvalues e - ^™)* of the pricing operator in 
the time changed model, where </>(A) is t he Laplac e expo nent of the Levy subordinator. Remarkably, this 
insight goes back to the seminal work of iBochnerl d 19491) introducing t he idea of subordination (see page 
370). It ha s been applied in probabi li ty the o ry (Albeverio and Riidigerl (120031): Chen and SonsJ d2005l)) and 
i n finance (lAlbanese and Kuznets"ov J 2004 Bovarchenko and Levendorskiil ( 20061 ): iMendoza-Arriaea et al 
(l2010l) : iMendoza-Arriaga and Linetskvl (|201 lh : [Li and Linetskv! (l201ll)). The idea t o use subordinated dif- 

> Ba 



fusions to build financial models goes back to iBarndorff-Nielsen and Levendorskiil (|2001l) , who considered 
NIG-like Feller processes for option pricing in the pseudo-differential operator framework. 

Surveys on the application of e igenfunction expansions to the valuation of European-style derivatives can 
be found in lLinetskvl (|2004l . 2008 ). where extensive bibliographies are given. Applications t o inter e st rat e 
models and the valuation of bonds without embedded options in p articular can be found in Lewis! 



Bovarchenko and Levendorskii 



1998 



2006). 



IDavvdov and Linetskvl (|2003l ): iGorovoi and Linetskvl (|2004 120071) : 
Applications of eigenfunction expansions to European-style derivatives i n models with jumps cons t ructe d 
by time changing diffusions with L e vy subordinators can be fou n d in lAlbanese and Kuznetsov (I2004I) . 
iBovarchenko and Levendorskiil ( 2006 ). iMendoza-Arriaga et al.l (120101). Mendoza-Arriaga and Linetskv 4 , 201 ll ). 
and lLi and Linetskvl (|201ll ). The reference IBovarchenko and Levendorskiil (|2006l ) is particularly relevant to 
our paper, as they also consider interest rate models based on subordinated diffusions. 

The rest of the paper is organized as follows. Section 2.1 describes the general framework for the 
application of eigenfunction expansions to short rate diffusion models. Section 2.2 describes short rate 
models with jumps constructed by time-changing diffusion models with a Levy subordinator. Section 3 
presents our eigenfunction expansion method for solving the dynamic programming backward induction for 
callable and putable bonds. Section 4 presents examples of eigenfunction expansions for CIR, subordinate 
CIR, Vasicek, subordinate Vasicek, the 3/2 model, and the subordinate 3/2 model. Section 5 presents 
numerical experiments demonstrating computational performance of the method on the test case considered 
in the previous callable bond literature. Appendix contains selected proofs. 



2 Short Rate Models 



2.1 Short Rate Diffusion Models 



Let {X t ,t > 0} be a one-dimensional, time-homogeneous regular (i.e. it reaches every point in (l,r) with 
positive probability) diffusion process on the interval ICR with (finite or infinite) endpoints I and r, 
— 00 < I < r < 00. An endpoint is unattainable if i t is a natural or an entrance b oundary and is attainable if 
it is an exit or a regular boundary (see pp. 14-15 of Borodin and Salminenl (2002) for Feller's classification of 



boundaries for one-dimensional diffusions). In this paper we assume that the endpoints are either unattain- 
able (and, thus, not included in the state space, so the interval / is open at such an endpoint) or regular 
and specified as instantaneously reflecting (the endpoint is included in the state space in that case, so the 
interval is closed at such an endpoint). We assume that the diffusion is conservative, that is P t (x,I) = 1 for 
each t > and x 6 /, where Pt(x, A) is the transition function from the initial state x to the Borel set AQ I 
in time t. Thus the process X has infinite lifetime. We assume that the volatility <r(x) of X is positive and 
continuous on the open interval (l,r) and the drift n(x) is continuous on (/,r). 

We further assume that the instantaneous interest rate (the short rate) r t at time t is a function of the 
state variable X t and is given by r t = r(X t ). We assume that r(x) is continuous on (l,r). The continuity 
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assumptions for a, fj, and r are not necessary, but simplify exposition in what follows. Consider the family 
of pricing operators (mathematically, Feynman-Kac (FK) operators) {VI, t > 0} defined by 



f(X t 



where M x is the expectation operator with respect to the probability measure ¥ x of the process X starting 
at x € I at time zero. Since in this paper we are interested in pricing, we always work with risk-neutral 
probabilities chosen by the market. The pricing operator V[ maps future payoff functions of the future state 
at time t into present value functions of the present state at time zero by discounting from time t to time 
zero and taking the expectation conditional on the current state at time zero. Under our assumptions, these 
operators form a strongly continuous semig roup on the Banach spa c e Cb( I) of bounded continuous functions 



on I with the supremum norm (see, e.g., iBorodin and Salminenl (120021) ) . The infinitesimal generator Q 



of this semigroup acts on twice continuously differentiablc functions on / with compact supports by the 
second-order differential operator (the Sturm- Liouville operator): 

G r f{*) = \o- 2 {x)f"{x) + n(x)f'(x) - r(x)f(x). 
Define s(x) and m(x) to be the scale and speed densities of the diffusion process X: 



s(x) = cxp 



My) 



dy } , m(x) 



x o 2 {y) J ' a 2 {x)s(x)' 
(l,r) is an arbitrary point (see Karlin and Tavlor ( 198lh . IBorodin and Salminen ( 2002 ) for 



where xq € 

discussions of the scale function and the speed measure of the one-dimensional diffusion). The infinitesimal 
generator can then be re-written in the formally self-adjoint form: 



G r f{x) 



1 



/'(*) 
m(x) \ s(x) 



r(x)f(x). 



Under our assumptions, the generator Q r and the FK semigroup {PI, t > 0} can be extended to a self-adjoint 
operator and the symmetric strongly-continuous semigroup in the Hilbert space L (7, m) of functions on / 
square-integrable with the speed measure m{dx) — m{x)dx and endowed with the inner product 

(f,g)= / f{x)g{x)m(x)dx 



and norm ||/|| = \/ (/, /)• Thus, the Spectral Theorem for self-adjoint operators in Hilbert spaces can now 
be applied to write down the spectral decomposition of the generator and th e semigro u p. T he spectral 
representatio n for one-dimensional diffusions goes back to the classical work of iMcKeanl (|1956l ). We refer 
the reader to iLinetskvl (|2004l [2008) for surveys of applications in finance. 

In this paper we further assume that the negative of the infinitesimal generator — Q r has a purely discrete 
spectrum in L 2 (I, m) bounded from below. Sufficient conditions for the purely discrete spectrum in terms of 
the behavior of the functions cr, /x and r near the boundaries / and r are given in lLinetskv (|2004ll2008h . When 
the spectrum of —Q r is purely discrete and bounded from below, the FK semigroup has the eigenfunction 
expansion of the form: 



VU{x) 



f n e Xnt ip n {x), f n = (f,(f n ), 



n=0 



(1) 



for any / G L 2 {I, m), where {A„}^ such that Ao < Ai < lim^oo A„ = oo, are the eigenvalues of 
— Q r and ip n are the corresponding eigenfunctions normalized so that ||<^n|| 2 = 1 (for future convenience we 
index the eigenvalues and eigenfunctions starting from n = rather than n = 1). The eigenfunctions form 
a complete orthonormal basis in L (I,m). We also assume that the eigenvalues satisfy the condition: 



n=0 



-A„t 



< OO 



(2) 
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for all t > 0, so that the FK semigroup is trace class (see Section 7.2 in IDaviesl ( 2007h ). We recall that the 
semigroup of a one-dimensional diffusion has a symmetric density pt{%, y) w ith respect to the speed measure 
m (dx) = m(x)dx that is a cont inuous function in t, x and y (see p. 149 in llto and McKeanl (|l974l ) or p. 13 
of Borodin and Salminen d2002h ). Hence, according to Theorem 7.2.5 in IDaviesl (]2007l ). the eigenf unctions 



tp n {%) are continuous functions with the global estimate |y> n (x)| 
density pt(x, y) has an eigenfunction expansion for alH > 



< 



e A ™*/ 2 y/p t (x,x) for all t > 0, and the 



Pt(x,y) 



OO 



-\ n t 



(3) 



that converges uniformly in x and y on compacts. This ensures that, in addition to the L 2 convergence, the 
eigenfunction expansion JT]) converges uniformly in x on compacts for all / <E L 2 (I, m) and t > 0. This follows 
from the Cauchy-Schwartz bound for the expansion coefficients < ||/||, the eigenfunction estimate, and 
the trace class condition ([2]). 

Since we are interested in bond pricing in this paper, we assume that the constant payoffs are in the 
Hilbert space L 2 (I,m), i.e. 1 e L 2 (I,m). This is equivalent to assuming that the speed measure m is a 
finite measure on /, m(I) < oo. In that case, the speed density can be normalized to one to be a probability 
density and, thus, serves as the steady state density of the underlying process X. Then the present value at 
time zero of a zero-coupon bond with unit face value and maturity t > when the underlying process is in 
state x, Xq — x, has the eigenfunction expansion given by: 



P(t,x)=E x 



-J*r(X u )du 



n=0 



-\„t 



(4) 



with the expansion coefficients p n — (l,ip n ). Under our assumptions, the expansion converges uniformly in 
x on compacts for all t > 0. 

Virtually all popular short rate models in the financial econ omics literature fi t into the gener al frame- 
work described above, including Cox-Inger soll-Ross (CIR) model ( Cox et al. . 19851) . Vasicek model (jVasicek . 
1977 ). the 3/2 model ( Ahn and Gaol 19991) . Black's model of in t erest rates as options dGorovoi and Linetskv . 
2004), and the quadratic model ( Beaglehole and Tennevl . [1992 ; Leippold and Wu . 2002 ). Note that we have 
not made the assumption that r(x) is non-negative to accommodate the Vasicek model. If we make that 
assumption, then Q r is positive semi-definite (so that — Q r is negative semi-definite), the FK semigroup 
{Vt)t>o is a contraction semigroup on L 2 (I,m), and Ao > 0. To accommodate the Vasicek model, we made 
a weaker assumption that the spectrum of —Q r is bounded from below, rather than non-negative. 



2.2 Short Rate Models with Jumps Constructed by Subordination 

A subordinator {T t ,t > 0} is a nondecreasing Levy process with the Laplace transform given by the Levy- 
Kchintchinc formula 

/•OO 

E [e- XTt ] = e -** (A) , 0(A) =7A+ / (l-e- A >(ds), A>0 



with the Laplace exponent 0(A), with nonnegative drift 7 > 0, and Levy measure i/(ds) satisfying the 
integrability condition J Q (s A l)z^(ds) < 00. For any set A C R bounded away from zero, jumps of 
sizes in A arrive according to a Poisson process with the arrival rate v{A). If v is a finite measure on 
(0,oo), the subordinator is a compound Poisson process plus drift at the rate 7. If i/(0,oo) = 00, the 
subordinator is a jump process with infinite activity and drift 7. If 7 = 0, it is a pure jump process. 
Examples of subordinators important in applications include compou nd Poisson processes with exponential 
or gamma distr ibuted jump s i zes, in verse Gaussian (IG) subordinators ( Barndorff-Nielserl 19981) . and gamma 



subordinators ( Madan et all Il998l) , with the latter two having infinite activity. These examples are special 



cases of subordinators with Levy measures of the form u(ds) — Cs p e ,,s ds with C > 0, rj > 0, and p < 1. 
The case with p e (0,1) are the tempered stable subordinators (the limiting cases with r\ = are stable 
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subordinators) . The special case with p = 1/2 is the IG subordinator. The limiting case with p — is the 
Gamma subordinator. The Laplace exponent is given by: 



0(A) 



7 A - CT(-p) [(A + riy - v p ] , p^o 



-/X + CHl + X/rj), 



p = Q 



(5) 



As an example, the Levy measure and Laplace exponent for an IG subordinator parameterized with [i and 
v i the mean and variance of an IG process at time one, t = 1, are given by: 



„(d*) = exp {-£.} ds, HX) = 7 A + £ 



2^A 



1 



(6) 



Further mat hematical details on sub ordinators can be found in Schilling et al. ( 2010( ) and on financial ap- 
plications in ICont and Tankovl (|2004r i . 

Since subordinators are non-negative, non-decreasing processes, they can be used as stochastic time 
ch anges to t ime change other processes. This procedure is known as Bochner's subordination and goes back 
to Bochner ( 19491 1955 ). In particular, time changing a Markov process with a subordinato r yields a nothe r 
Markov process whose s emig roup and infinites imal generator are giv en by Phillips' theorem ([Phillips! (|1952[ ). 
Theorem 32.1 in Sato! ( 19991) . Chapter 12 in Schilling et al.l (l201(t)') : For recent financial a pplica tions see 
Mendoza-Arriaga et al. ( 2010l ). Mendoza-Arriaga and Linetskvl ((201 lh . and iLi and Linetskv ( 2011 ). 

In particular, we can construct new short rate models with jumps from diffusion short rate models as 
follows. Let X be an ergodic one-dimensional diffusion on / with volatility o~{x) and drift fi(x) and with the 
stationary density given by the (normalized) speed density m as described in the previous section. Let r{x) 
be the function defining the diffusion short rate model as in section 2.1 and {VI, t > 0} the corresponding 
FK semigroup. Then a new short rate model with jumps is obtained by subordinating the original diffusion 
short rate model with respect to a given subordinator T with the Laplace exponent cj)(X) by defining a new 
semigroup (the superscript <f> signifies that the subordination is performed with respect to the subordinator 
with the Laplace exponent 4>): 



V r t '+f(x) :=E X \e-ti^)*u f{ y t) 



t > 0, 



(7) 



where Y t is a new Markov jump-diffusion process on I with the infinitesimal gener ator acting on twice- 



differe ntiab le functions with compact suppo r ts as an integro-diffcrcntial operator (see lMendoza-Arriaga et al 
(|201dh and iMendoza- Arriaga and Linetskvl (|201lh ): 



/(*) = \l°*V)f"V) + li*(x)f'(x) + ^ (f(y) - f{x) - l { \ y - x \< l} {y - x)^(x)j ^(x, y) dy, 



where the drift with respect to the truncation function l/i^-x^i} is 



^{x) = !H{x) + / / {y - x)p s {x,y)Ay \v{As), 

J(0 : oo) \J {ye(l,r):\y-x\<l} ) 

and TT^(x,y) is the symmetric state-dependent Levy density: 

7r*(x,j/) = / p s {x,y)v{ds), 

where p s (x,y) is the density of the original FK semigroup (Vl)t>o of the pure diffusion short rate model, 
and 7 and v are the drift and the Levy measure of the subordinator. When 7 > 0, Y is a jump-diffusion. 
When 7 = 0, Y is a pure jump process. Y has the same steady state density m(x) as the original diffusion 



G 



X. The short rate in Q is the function of the state variable Y t so that r t — r^(Y t ) with the function r^(x) 
given by: 



r*(x) = ~/r(x) + / (1- P(s,x))v(ds), 

J(0,oo) 

where P(s, x) is the price ((4J of the s-maturity zero-coupon bond at time zero when the state of the underlying 
diffusion Xo = x. The form of the generator of Y and the function follow from Phillips' theorem that 
characterizes the subordinate semigroup and its infinitesimal generator. Here we subordinate the pricing 
(FK) semigroup of the diffusion process X with the discount rate r t — r{X t ) with respect to a given 
subordinator to obtain a new semigroup interpreted as the pricing (FK) semigroup of a new short rate 
model r t — r^(Y t ) driven by the state variable Y t following a Markov process with jumps. 

Remark 2.1. Mathematically, the subordination of the FK semigroup can be interpreted as follows. First 
formulate the FK semigroup V r of the original conservative diffusion X with the discount rate r(x) as 
the transition semigro u p of t he diffusion X with killing at the rate r(x) (cf. Section II. 4 on pp. 27-28 in 
Borodin and Salminenl ( 2002 ) for the connection between discounting and killing). Then construct a new 



process Xf := X% by time changing X with the subordinator T. Use Phillips' theorem to write down its 
infinitesimal generator and, thus, its local characteristics (diffusion, drift, state-dependent Levy measure, 
and state-dependent killing rate r*(x)). Then formulate the transition semigroup of the process X^ as 
the FK semigroup of a conservative process Y with the generator given above and with the discount rate 
r^jx) given above. The form ulation of the application of Phillips' theorem to this situation is given in 
Mendoza-Arriaga et al. ( 201dh and we do not repeat it here to save space. 



The subordinate FK semigrou p {P^ , t > 0} is also a symmetric strongly continuous semigroup of 
operators on L 2 (I, m) ( Chenl . l2005T ). and, under the assumptions we have made about the semigroup {VI , t > 



0} in the previous section, it possesses an eigenfunction expansion in the same eigenfunctions (p n (x) with A„ 
in Eq.|T]) replaced with A* := <j>(\ n ), where 4>(X) is the Laplace exponent of the subordinator: 

oo 

P?*f(x) = J2 /ne-^VnW, fn = (/, ¥>«), (8) 
n=0 

for any / € L 2 (J, m) and t > 0. For mathemati cal details on subordin ation of semigroups of operators and 
Ma rkov processes see the excellen t exposition inlSchilling et al] (120101) . an d for recent financial a p plicat ions 
see iMendoza-Arriaga et al. ( 2010l ). iMendoza-Arriaga and Linetskvl (| 2 1 lh . and Li and Linetskvl ( 201ll ). If 



we further assume that the Laplace exponent <f> of the subordinator is such that it satisfies the condition 

( A ")' < oo (9) 



71=0 

for all t > 0, then the semigroup {"P t r '^, t > 0} is also trace class. If we further assume that the eigenfunctions 
of the original pure diffusion FK semigroup have a bound independent of n on each compact interval K = 
[a,b] C {l,r), i.e, |<^„(x)| < Ck for all n, where the constants C'k may depend on the interval K but are 
independent of n, then these assumptions ensures that, in addition to the L 2 convergence, the eigenfunction 
expansion of the subordinate semigroup ([5]) converges uniformly in x on compacts for all / € L 2 {I,m) and 
t > and that the subordinate semigroup {Vl'^, t > 0} also has a continuous density with respect to m(x)dx 
with the eigenfunction expansion for all t > 

oo 

pt(x,y) = J2e-* (Xn)t Mx)<Pn(y) (10) 

n=0 

uniformly convergent on compacts in x and y. 

For the present value at time zero of a zero-coupon bond with unit face value and maturity t > when 
the underlying state variable has initial value of Yq — x, we then obtain the eigenfunction expansion given 
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by: 



p(t,x) = E x e -/o^)d« = y; Pn e-^«)Vn(x) (n) 



/ , j 

n=Q 



with the expansion coefficients p n — (l,(p n ). Under the assumption Q on the growth cj>(\) and the bound 
on eigenfunctions independent of n, the expansion converges uniformly in x on compacts for all t > 0. 

Using the subordination approach, we can extend all the diffusion short rate models popular in financial 
economics to jump-diffusion and pure jump models, in particular constructing subordinate CIR (SubCIR), 
subordinate Vasicek (SubVasicek), etc. Subordinate models allow for jumps in the interest rate dynam- 
ics. Moreover, if the diffusi on process is mean-rev erting, the subordinate process will have jumps that 
are also mean-reverting (see iLi and Linetskv ( 2011 ) for the proof in the subordinate Ornstein-Uhlenbeck 



context). While adding jumps improves the model's realism and flexibility, remarkably, the analytical and 
computational framework remains entirely unchanged, as the only modification required in the eigenfunction 
expansion is the replacement of A„ in Eq.((T]) with A^ = </>(A n ) in Eq.©. 

Remark 2.2. Matching the initial yield curve. The time-homogeneous short rate models discussed in sec- 
tions 2.1 and 2.2 c a n be extended to match any initial term structure of interest rates as proposed by 
Brigo and Mercuriol ( 2001 ) and commonly done in fixed income market practice by adding a deterministic 



function of time to the short rate process in the extended diffusion and the subordinated diffusion models, 
respectively. The function can then be chosen so that the initial zero-coupon bond prices of all maturities in 
the extended model match the zero-coupon prices consistent with the given initial term structure of interest 
rates (given yield curve) . The callable and putable bond pricing developed in this paper can then be imme- 
diately extended to this class of models. To simplify notation, we do not explicitly consider this extension 
in what follows and assume the short rate model is time homogeneous. 



3 The Eigenfunction Expansion Method for Callable and Putable 
Bonds 

The call option allows the bond issuer to buy the bond back from the bond holder (call the bond) for pre- 
specified call prices at some pre-specified times prior to maturity. This allows the bond issuer to refinance 
the bond if interest rates decline. The put option allows the bond holder to sell (put) the bond back to 
the bond issuer for pre-specified put prices at some pre-specified times prior to maturity. This allows the 
bond holder to re-invest the proceeds into a bond with higher coupon if interest rates rise. We assume that 
the bond principal is equal to one dollar and the bond pays coupons of C dollars on the coupon dates ti, 
i = 1, k. Let the bond issue date and maturity date be to = and tk = T, respectively. After some initial 
protection period from to until t^* , the call and put options can be exercised at the subsequent coupon dates 
prior to maturity tk ■ We also assume that there are notice periods of lengths 8 so that the 
option exercise decision to exercise at the coupon date U has to be made at an earlier time r» = U — 8 so that 
the adequate advance notice of duration 8 can be given to the other party. Denote the call and put prices at 
times ti as Kf and Kf, respectively. It is assumed that Kf > Kf. This sets up an optimal stopping game 
in discrete time with finite horizon where one player (the bond holder) chooses a stopping time to maximize 
the bond value and the other player (the bond issuer) chooses a stopping time to minimize the bond value. 

We assume that the short rate is r t = r^(Y t ), where Y t is a subordinate diffusion, as discussed in section 2. 
We note that the pure diffusion model can be viewed a special case with the trivial time change % — t with 
4>{X) = A, 7 = 1 and v = 0. We assume that all the assumptions made in section 2 are in force. Then the 
pricing operator has the eigenfunction expansion ([lip. To simplify notation, we drop the superscripts 
r and <f> and simply write Vt for the pricing operator in what follows. 

Let V(t, x) be the value of the bond at time t € [0, T] when the underlying state is x. Since the decisions 
to exercise the call and put options are made at times n = ti — 8 prior to the coupon dates, the present 
values at Tj of the call and put prices at time ti must be compared to the holding (continuation) values of 
the bond at time Tj. The discounted value of the call price at time Ti is KfP(8,x), where P(8,x) is the 
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value of the zero-coupon bond with time to maturity equal to the notice period S and unit face value when 
the state variable is in the state x at time r». The expected discounted value of the put price at time Tj is 
KfP(S, x). At maturity the bond value is equal to its principal plus the last coupon, V(tk, x) = 1 + C. Let 
V k (x) :— V(tk,x) and, for i < k — 1, V l (x) := V(ri,x) denote the bond's value at time Tj and C l (x) denote 
the bond's holding (continuation) value at time Tj ex-coupon at time ti, i — k* , k — 1 (the assumption is 
that the next coupon C is always paid at time U, whether or not decisions are made to exercise any of the 
two options at time Tj). Let V"°(a;) := V(0,x) denote the value of the bond at the time of issue. 

Assuming both players behave rationally and maximize the value of their assets and minimize the value 
of their liabilities, in state x the bond issuer will exercise the call option at time r» if K?P(6,x) < C l (x), 
the bond holder will exercise the put option at time Tj if KfP(6,x) > C l (x), and there will be no exercise 
of either option at time r, if KfP(S,x) < C l (x) < KfP(S,x). By the assumption that K? > Kf the 
simultaneous exercise of both options is never optimal. Further set hi := r»+x — Ti, i — k*,...,k — 2, and 
hk-i := tk — Tk-i- Then the value of the bond wit h call and put opt i ons sa tisfies the following Bellman's 



dynamic programming backward induction (see also iBen- Ameur et al.l ((2007)): 



V k (x) = (l + C), 

C l (x) =V hi V i+1 (x), k*<i<k-l, 
V\x) = max {K?P (5, x), min {K?P (6, x),C l (x)}} +CP(S,x), k* < i < k - 1, 

fe*-i 

V°(x) = V Th , V k ' (x)+CJ2 x )- 



(12) 

(13) 
(14) 

(15) 



Assuming that for each i = k*, k — 1 each of the two equations 

KfP(5, x) = C l (x) , KfP(5, x) = C l [x) 



(16) 



has at most one solution in / denoted by x\ and xf, respectively, setting x\ := I if the first equation has 
no solution in /, setting x\ := r if the second equation has no solution in /, and noting that x\ < x\, the 
backward induction (|14[) can be re-written in the form: 



V\x) = K?P(5, x)l {x<x c } + KfP(S, x)l {x><] + C l (x)l {x ^ x ^ } + CP{5, x), 



(17) 



for all k* < i < k — 1 and where li.y denotes the indicator function, and l{ x <x c } = if x\ — I and 

l{x>x p } = if x\ = r by convention. 

This backward induction can be solved by a variety of computational methods in the literature, as 
discussed in the introduction and in section 5, based on the different methods to approximate the pricing 
operator Vt appearing in Eqs. p^l) and ([T5"j) . In this paper we follow the approach based on representing the 
pricing operator Vt by its eigenfunction expansion. Our main result is the following theorem that summarizes 
our eigenfunction expansion method for the valuation of callable and putable bonds. 

Theorem 3.1. Suppose that m is a finite measure on I , m(I) < oo. 

(i) The value function V (x) and the value functions V i (x) and the continuation value functions C r {x) 
are in L 2 (I,m) for all i = k* , ... , k — 1 . 

(ii) The continuation value functions have the following eigenfunction expansions: 



oo 

C\x) = c^e-^cfnix), i = k*, k - 1. 



The value function at the time of the bond issue has the following eigenfunction expansion: 



oo k* — 1 

V°(x) = cfe- A ^>„(x) + C ]T P(U, 



x). 



(18) 



(19) 
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(Hi) The eigenfunction expansions (|18[) and Q19p converge uniformly in x on compacts in I. 
(iv) If each of the two equations in (|16[) has at most one solution in I, then the eigenfunction expansion 
coefficients in (|18p and (|19[) satisfy the following backward recursion: 

c k n = (l + C) Pn , n = 0,1,2,..., 

where p n = (1, cp n ) are the expansion coefficients of the unit payoff appearing in the eigenfunction expansion 
of the zero-coupon bond (fTTj) , and for each i = k*, k — 1, 

oo 

4 = i^„(Z, ij) + c^e-^TT^nixt, x?) + Kfp n (xl r) + Cp n e-^ 5 , (20) 

m=0 

/or n = 0,1,2,..., where x\ and x\ are as previously defined, and for I < x < y < r we introduced the 
following notation 

rv 

K m ,n{x,y) ■= 0-(x,y)<Pm,Vn) = / Pm {z)fn (z)m(z) dz, (21) 

J X 

Pn(x,y) := (l {x , y )P(6),(p n ) = ( P(5,z)cp n (z)m(z)dz, (22) 

J x 

where l( x ,y) = l( x .y)( z ) *- s ^ e indicator function of the interval (x,y) and P(t) = P(t,x) is the value function 
of the zero-coupon bond with time to maturity t > and unit face value when the underlying state variable 
is x. 

Proof. See the Appendix. □ 



Theorem 3.1 reduces the solution of the backward induction for callable and putable bonds to the recursion 
for the expansion coefficients (|2"0"1) . together with finding the roots of equations (| 16[) by a numerical root 
finding algorithm, such as bisection. The continuation value C l {x) on the right hand side of the equations 
(|16p is given by the eigenfunction expansion (|18|) with the coefficients determined on the previous step of 
the recursion. The eigenfunction expansion can be truncated at a finite level, and the truncation error is 
uniformly controlled due to the uniform convergence of the expansions when condition ([SJ is satisfied. 

Remark 3.1. When the state process is an affine diffusion, such as CIR or Vasicek, the zero-coupon bond 
value function has the exponential-affine form in the state variable 

P(t,x) = A(t)e- B ^ x , (23) 

and the integral in Ea. (|22|) can be written as 

rv rv 
p n (x,y)= P(t,z)<p n (z)m(z)dz = A(t) e- B ^ z <p n (z)m(z)dz (24) 

J X J X 

and in some cases can be calculated in closed form. Generally, it can be calculated using the eigenfunction 
expansion for the zero-coupon bond value function: 

oo 

(l( x , v )P(t),<p n ) = Pme~ HXm)t K m ,n(x, y), (25) 
where Tr m ,n(x, y) are defined in Eq. (j2"T|) . 



Remark 3.2. If the bond has only the call option and no put option, then the game reduces to the optimal 
stopping problem for the bond issuer. In that case in Ea. (|2"0"|) the term with Kfp n (x l p ,r) is absent and 
Tr m .n(Xi , x?) is replaced with 7T TOin (x?, r) for all i. Similarly, if the bond has only the put option and no call 
option, the game reduces to the optimal stopping problem for the bond holder and in Ea. (j20"]) the term with 
Kfp n (l,x l c ) is absent and ■n m ^ n (x1,x p i ) is replaced with ir m<n (l, x?) for all i. 
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The condition that each of the two equations in (|16p has at most one solution in I generally needs 
to be checked case by case. For CIR, Vasicek, and 3/2 models considered in section 4 the condition can 
be verified by the fo l lowing proposition. For CIR and Vasicek, (|26p below follows from Theorem 1.1 in 



Ikeda and Watanabe (|l977h while holds for the 3/2 model since the diffusion process X(t) in the 3/2 



model can be written as X(t) = 1/Y(t), where Y(t) is a CIR process satisfying the Feller condition. 

Proposition 3.2. Suppose that the discount rate r(x) is a non- decreasing function. Let (Q, J-, P) be a 
complete probability space with right continuous increasing family (J t ) ( >o of sub a-fields of J- each containing 
P-null sets and let B t be a one- dimensional J- t -Brownian motion. Let o~(x) and ^(x) be continuous. Let 
Xi(t) and X2{t) be processes started at different initial states such that X\{Q) < XzifS) and 

Xi{t) = Xi(0) + f a(Xi(s)) dB(s) + f n(Xi(s)) ds, i = 1,2. 
Jo Jo 

If we have that 

F[Xt(t) < X 2 (t) for allt>0] = l, (26) 



then each of the two equations in (|16[) has at most one solution in I . 

Proof. See the Appendix. □ 



Remark 3.3. We emphasize that uniqueness of roots of (|16[) is not a requirement for our method to work. 
In fact, it is one of the strengths of our approach that it can handle just as easily more general cases 
with multiple break-even points and, hence, early exercise regions that are not necessarily one-sided and, 
in general, can be unions of multiple intervals. Proving the uniqueness of the break-even point provides a 
convenience for numerical implementation, as we can stop after finding a single root. In general, without 
the proof of uniqueness, a more thorough numerical investigation of the functions is required in each case 
to either establish uniqueness or determine multiple roots. While we have been able to prove in Proposition 
3.2 uniqueness for pure diffusion short rate models, we have been unable to extend the proof to the case of 
subordinate diffusions, as it is based on classical SDE comparison results that to the best of our knowledge are 
not available in general for subordinate diffusions. Nevertheless, in our extensive numerical experimentation 
in all cases of subordinate diffusions we have considered, we have observed similar behavior of functions in 
(|TB)) that lead to unique solutions. We thus conjecture that uniqueness also holds for subordinate diffusions, 
perhaps subject to some condition on the subordinator. 



4 Examples 

4.1 CIR and SubCIR 

In the Cox-Ing ersoll-Ross (CIR) model (|Cox et all Il985l) . the short rate follows the CIR diffusion with 



drift (i(x) — k(9 — x) and volatility cr(x) = o-\fx where k > 0, 9 > 0, and a > are the rate of mean 
reversion, the long run mean and volatility, respectively. In this case r(x) = x. When Feller's condition 
2nd I 'a 1 > 1 is satisfied, the origin is an unattainable entrance boundary and infinity is an unattainable natural 
boundary. In this case / = (0, oo). When Feller's condition is not satisfied, the origin is an attainable regular 
boundary and is specified as instantaneously reflecting. In this case / = [0, oo). The CIR speed density reads 
m(x) = -^2X b ~ 1 e~~^ . 

The celebrated CIR zero-coupon bond pricing formula is: 

P(t,x) = A(t)e- B ^ x 7 (27) 

where 

*m-( -. , 2 r':r"" V. m- 2(e7 '- 1) 



(7 + /c)(eT*- l) + 2 7 y ' w ( 7 + K)(e'>' t - 1) + 2 7 ' 
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The eigenfunction expansion (j4]) of the CIR zero-coupon bond pricing function is given in lDavvdov and Linetskv 



(2003). In this case the eigenfunctions, eigenvalues, and the expansion coefficients for the unit payoff arc: 



\ 9 / ' 

V a ) 



N n = 



2j\ 

2T(b + n) \a 2 J 



fa/2 



(1, <Pr. 



2N n T(b + n) ( a 



7 + K 



k — 7 
K + 7 



where L„ (a;) are the generalized Laguerre polynomials and b and 7 as defined in 

The CIR eigenfunctions are continuous and have a bound independent of n on each compact interval 
K = [a, b] C (I, r ), i.e. |(p n (x)| < Ck fo r all n, where the constant Ck is independent of n, since by inequality 
(27a) on p. 53 of lNikifo rov and Uvarovl (|l988h . for any compact interval in J, the CIR eigenfunctions satisfy 
the bound |y n (a;)| < Crt~ 1//4 , where the constant C is independent of n (but depends on the interval). 
The quantities ([2"T]) and (j2"2")) in the CIR model can be calculated as follows: 



2^J 



7 

b-l 



,(6-1) 



U 2 



,(6-1) 



(2jX 



(29) 
(30) 



where we introduced the following notation: 



B{S)a 2 



K + 7 



2 7 27 ' 

a^l(x)= I L^(y)L^(y)e-yy a dy, bi a \s,x) 



In the calculation of ([30)1 we used the explicit expression for the CIR zero-coupon bond pricing function ([27|) 
as in Eq. (|24p . rather than its eigenfunction expansion. 

The quantities a[ a i(i) and b„ (s,x) with a = b — 1 > —1 can be efficiently computed via the following 
recursion. 

Proposition 4.1. Suppose that a > —1. The coefficients al^m(x) are computed as follows for all x > 0. 
For n > 1, m> 1, m ^ n, 



For n > 1 , 



For m = n. 



aH(x) = l -e~*x^L^\x). 



a$ = 7 (a + l,x), a%l(x) = \ [L£\x)L<g L \x)e-*x a+1 + a^^x) 

where 7(0: + 1, x) = J? e~ v y a dy is the lower incomplete gamma function. 
The coefficients (s, x) are computed recursively as follows for all x > 0. 



n > 1, 



(31) 
(32) 

(33) 



= -^(a + Mz), b^(s,x) = i e — ^I^W 



^e^W, n>l, 



(34) 
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Proof. See the Appendix. □ 

The Laguerre polynomials of degree a satisfy the following classical recursion (<Lebedevl . [l965l Eg. 4.18.1) 

a-l-x\^ a) ( a-l\ (o) 



and 4 Q) (a;) = 1, L[ a \x) 



—x + a + 1. Then for any TV the quantities {ai a „(a;),0 < n < N} can be 
efficiently computed recursively in the following order: 



• a. 



(a+N) 
0,0 



(x) 



• a. 



0.0 



'■o 



• a O,o( 2 '): a Ll( x )i a 7V-l,JV-l( x )' a 7V,Jv( 2 ') 

The computation of al"i(a;) with n ^ m can be done directly using (1311) and (1321) and the recursion for the 
Laguerre polynomials. The quantities {bl?\x),0 < n < N} can be efficiently computed recursively in the 
order: 



i i a N-l,N-l\ 



(«) 



» 



• b, 



(a+N) 




(X) 



d<*+l) 



(x),b[ a+1 \x),...,b^>(x) 



. b^(x),b i r\x), ■■■,b { °Ux),b i £\x) 

For the SubCIR model the explicit bond pricing formula similar to ((27)) is not available, and we use the 
eigenfunction expansion (fTTj) of the SubCIR zero-coupon bond price as in Eq. (|2"5|) . The expression (|3"0)) is 
then replaced with: 



Pn(x,y) ^^Pr, 



-0(A m )«5 



m=0 



'2l 



6-1 



7 



,(6-1) 



V a 2 



,(6-1) 



/27X 



(35) 



In Theorem 3.1 A„ are now the eigenvalues of the SubCIR model related by </>(A n ) to the eigenvalues of the 
CIR diffusion model. 

Remark 4.1. In the limiting case x = 00 we have 

7 f^l 



N n N m V a 



T(a + n + l)(s - l) f 

nloa+n+l 



due to the orthogonality of Laguerre polynomials and the integral identity (jGradshtevn and Rvzhild . 12007 , 
p.809) 

where a > — 1, s > 0, n = 0, 1, 2, .... Using these coefficients, the recursion for the expansion coefficients in 
Theorem 3.1 simplifies in the case of callable bonds with no put option (in that case Xi — 00 and there is no 
term with Kf in Eq.(20)). 



13 



4.2 Vasicek and SubVasicek 

In Vasicek model (jVasicek , 19771) . the short rate follows the OU diffusion with drift fi(x) — n{9 — x) with 
K > and 9 > and constant volatility a > 0. In this case / = M and r(x) = x, both boundaries at 
plus and minus infinity are unattainable natural boundaries, and the process can get negative. However, 
when 9 and the initial state xq are sufficiently above zero and k > is sufficiently large, the probability of 
the rate falling below zero is relatively small due to mean reversion pulling the process back towards the 
positive long run mean as it approaches zero from above. The Vasicek speed density is a Gaussian density 

k(9-x) 2 

m(x) = -jze . 

The celebrated Vasicek zero-coupon bond pricing formula has the same exponential afhne form as the 
CIR (EH) with 



B(t) 



(1 - e- Kt ) , A(t) = exp 



K z 4k 



The eigenfunction expansion (j4|) of the Vasicek zero-coupon bond pricing function is given in lGorovoi and Linetskv 



(|2004l ). In this case the eigenfunctions, eigenvalues, and the expansion coefficients for the unit payoff are: 

2 



A, 



a" 
2^? 




k a 



n2 n+1 nV 



Vn {x)=N n e'<-^H n {i + a), £ : =^(a:-0), a —, N n -- 

a k ' 

2 Fir _o£ 
Pn = -\ -N n a n e * , 

(7 y K 

where H n (x) are Hermite polynomials. 

The Vasicek eigenfunctions are continuous and have a bound independent of n on each compact interval 
K = [a, b) C (I, r), i.e. < Ck for all n , where the constant Ck is independent of n, since by inequality 

(28a) on p. 53 of lNikiforov and Uvarov ( 1988 ). for any compact interval in /, the Vasicek eigenfunctions satisfy 
the bound |<^„(a;)| < Cn -1 / 4 , where the constant C is independent of n. 

The quantities (l2~Tj) and (|2"2"|) in the Vasicek model can be calculated as follows: 



. 2A{5)N n . 
Pn{x,y) = -;=- e 



- B (&)(e-^L) 



'».., ! s, — (y - 0) + a 



where we introduced the following notation: 



B(S)a 



— (x -9) + a 
a 



I),, ( s, (x — 9) + a 



(36) 
(37) 



e v H n (y)H m (y)dy, b n (s,x) := 



H n (y) dy. 



In the calculation of p7p we used the explicit expression for the Vasicek zero-coupon bond pricing function 
as in (|24j) . rather than its eigenfunction expansion. 

The quantities a„ j?n and b n can be compute d efficient l y. Th e coefficients of the Hermite polynomial H n (x) 
can be computed from the recursive equation ( Lebedev . 19651 Eq. 4.10.1) 



H n (x) = 2xH n _ 1 (x) - 2(n - l)H n _ 2 (x), n > 2, H (x) = 1, H x {x) = 2x. 
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Proposition 4.2. The quantities a rn , n (x) can be computed as follows. For m 7^ n 

H n (x)H rn+1 (x) - H m {x)H n+1 (x) _ x 2 

a n ,m{x) = —. r e . (38) 

2{m — n) 

ct n ,n{x) can be computed recursively as follows: 

a ,o(x) = y/TT$(V2x), a n>n {x) — -H n ^ 1 {x)H n (x)e~ x + 2na„_ 1)n _i(a;), n > 1, (39) 

where <&(x) is the standard normal cumulative distribution function. 
The quantities b n (s,x) can be computed as follows. 

&o(*,z) = ^e£^(ErfQ(2i-«))+lJ, (40) 

where Erf (x) is the error function, and 

b n (s,x) = - e sx ~ x -ff„_i(x) + sb n -x(s,x), n > 1. (41) 

Proof. See the Appendix. □ 

For the SubVasicek model the explicit bond pricing formula similar to (|23p is not available, and we use 
the eigenfunction expansion (jllj) of the SubVasicek zero-coupon bond price instead. The expression (|37p is 
then replaced with: 



o n (x,y) Pi 



ro=Q 



1 rr 



a n , m I — (y - 0) + a ) - a n , TO ( — (x - 6) + a 



(42) 



In the recursion ([2U|) A n are the eigenvalues of the SubVasicek model related by <fi(\ n ) to the eigenvalues of 
the Vasicek diffusion model. 

Remark 4.2. In the limiting case x = 00 we have 

a m , n (oo) = ^ (5 m ^„, 6„(s, 00) = e s2/4 \/7i : (-s)" 

ziv n iV m 



due to the orthogonality of Hermite polynomials and the integral identity ( Gradshtevn and Rvzhikl . 2007 
Eq. 7.374.6, p.803) 

/•DO 

-^ 2 i/„(x) &x = ^(~2y) n . 



e 
—00 



These coefficients can be used to evaluate callable bonds without the put option, similar to remark 4.1 for 
the CIR. 



4.3 The 3/2 and Sub-3/2 Model 

In this model, the short rate process is a diffusion on (0,oo) with infinitesimal parameters o~{x) = o~x 3 / 2 , 
n( x) = k(9 — x)x, r (x) = x, where k, 9, and a are pos itive constant param eters. This process was proposed 
bv lCox et al. I (Il985h as a model for the inflation rate, n and Gaol (Il999h propose this process as a model 



for the short rate and show that this model is empirically more plausible than the square-root model. 
Let a = At + 1, = m = J + \) 2 + ^i- The speed density for this model is m(x) 



2a 1 e x. The eigenfunction expansion (U) of the zero-coupon bond pricing bond is given in iLinetskv 



(|2004l ). In this case the eigenfunctions, eigenvalues, and the expansion coefficients for the unit payoff arc: 

A n = n9(n + m — a + 1/2), 
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<p n {x) = N n x a - m - l ' 2 L^ (tj , N n = ^, 



2f{2m+l n \ 



O- 2 



2T(2m + n+ 1)' 

_ Aw fl -a-m-* T(g + m + l/2)r(m + n - a + 1/2) 
Pn-ga^P - nW(m-a + 1/2) 

The 3/2 eigenfunctions are continuous on I and have a bound independent of n on each compact interval 
K = [a, b) C (I, r ), i.e. |</? ra (x)| < Ck for all n. where the constant Ck is independent of n, since by inequality 
(27a) on p. 53 of iNikiforov and Uvarovl ( 19881) . for any compact interval in I, the 3/2 eigenfunctions satisfy 
the bound |y>„(x)| < Cn -1 ' 4 , where the constant C is independent of n. 

The quantities (l2~Tj) and (|22p in the 3/2 and Sub-3/2 model can be calculated as follows: 



Pn(x,y) = ^Pke 



2N n N k 

CT 2^2m+l 
-X k 5 ™ n N k 



(2m) ( P 



- I - a. 

X 



(2m) ( P 
k 



k=0 



(2m) 
l n,k 



(2m) / P 



a 2 P 2m+1 

where we introduced the following notation: 



n . k 



(43) 
(44) 



For the Sub-3/2 model, A„ are the eigenvalues of the Sub-3/2 model related by </>(A ra ) to the eigenvalues of 
the 3/2 diffusion model. 

The expression for a n ^k(x) is same as the one for the CIR model. Hence, the quantities a nj k(x) can be 
computed using the method given in the CIR section. 



5 Computational Results 



This section shows our computational results for CIR, Vasicek, SubCIR, and SubVasic ek models. We consider 



the ca llable bond exa mple that has been ext e nsively used in the l itera ture star ting fromlButtler and Waldvogel 
(|l996l ) and including kTHalluin et ail (|200l[ ). lBen-Ameur et ail ([20071 ). and lde Frutol ([20081 ). as the test case 
to compare computational performance of a number of computational approaches to the callable bond val- 
uation. The callable bon d was issued by Swiss Confede ration in 1987 with maturity in 2012. At the time 
of valuation considered in Biittler and Waldvoeell (|l996l) and the subsequent papers, the remaining time to 
maturity of the bond was t k = 20.172 years with k — 21 remaining annual coupons of 4.25% per annum. 
The notice period is 2 months, 5 = 0.1666. The protection period is tk* = 10.172 with k* = 11. There are 
ten early exercise dates, t\x — 10.172, t\% — 11.172,..., t 2 o = 19.172. The call prices corresponding to these 
dates are given in Table Q] The bond did not inc lude a put option. We use the values of the parameters n, 
a, 9 for the CIR and Vasicek models estimated in lButtler and Waldvoeell ( 19961 ) and used in the subsequent 
papers in the literature. They are given in Table [5] For the SubCIR and SubVasicek models, we used the 
same parameter values for the underlying CIR and Vasicek diffusions, while specifying the subordinator to 
be the inverse Gaussian (IG) subordinator with drift (IG Levy measure given in ([5])). The IG parameter v 
was set to 1. The subordinator drift 7 and the IG parameter \x were chosen so that E[T t ] = t to normalize 
the time change. For the jump-diffusion case, we used 7 = 0.5 and /1 = 0.5. For the pure jump case, we used 
7 = and fi = l. 

In the process of finding the break-even point x\ at each step of the recursion, the infinite series for the 
continuation value given in ([T%|) needs to be truncated at some finite level. At time to, the series in (TTH]) 
also needs to be truncated. In the recursion formula for the expansion coefficients in (f2"0"f , only the previous 
expansion coefficients that were calculated at the previous step are used in the eigenfunction expansion. For 
the subordinated models, we also compute the zero-coupon bonds by the eigenfunction expansions that are 
also truncated at some finite level. We used an adaptive truncation strategy that truncated the expansion 
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after a user-specified relative error tolerance e was reached in each instance of the series evaluation by 
comparing with e the ratio of the next term and the sum of the next two terms relative to the sum of all the 
previous terms. 

In order to find the break-even point at each decision point in time, the bisection method was used. 
The break-even short rates are shown in table [3] (JD stands for jump-diffusion, and PJ for pure jump). For 
Vasicek and CIR models, r(x) — x, so the short rate is equal to the state value. For subordinated models, 
the short rate is given by the function r^(x) of the state variable. For the CIR model, we start by checking 
the boundary at zero to see if the expected discounted value of the strike is greater than the continuation 
value, in which case there is no non-negative break-even point at that decision time instance. Otherwise, 
there is a unique non-negative break-even point for the CIR model. For the Vasicek model, there always is 
a unique break-even point. The break-even point was found by the bisection method until the length of the 
search interval became less than 10~ 7 . 

Table |4] shows computational results for all the models considered in this section with the initial short rate 
r = 0.05. The first column indicates the absolute pricing error in pricing the callable bond. During the process 
of finding the break-even point, the truncation level is determined for each evaluation of the continuation 
value. The second column gives the average truncation level N at each decision point at time T20, rig, rn, 
and at to (the average of truncations levels as determined by our adaptive truncation algorithm in evaluating 
the expansion of the continuation value needed for each step of the bisection algorithm). The third column 
shows the maximum truncation level N at those times. The fourth column shows the CPU time of our 
algorithm implemented in C using the GNU Scientific Library (GSL) and compiled with gcc and executed 
on a 2.4 GHz Intel Core i3 370M processor. The CPU time includes the time taken for any precomputations 
required. For CIR and Vasicek models, the CPU time to price the callable bond to approximately five correct 
decim al points (convergence of 10~ 5 in the tables) was about one millisecond. F or comparison, de Frutosl 



(2008) reported CPU times of 0.75 seconds using Matlab on a 3GHz processor, while Ben-Ameur et al. (2007) 



reported CPU times of 2 to 3 seconds using C on a 2.0 GHz Pentium 4 processor. Thus, the eigenfunction 
expansion approach is approximately three orders of magnitude faster in this instance of pricing the callable 
bond. 

For the subordinated models with jumps, jump-diffusion models required slightly longer CPU times than 
pure diffusion models, while pure jump models required slightly longer times than jump-diffusions. This is 
due to the replacement of the diffusion eigenvalues A„ with eigenvalues </>(A n ) of subordinated processes that 
slows down the eigenvalue growth and, hence, requires more terms in the expansions, as evidenced in Table 4. 
Still, the algorithm reached the pricing error of under 10 -5 in 2.2 and 2.5 milliseconds under pure jump CIR 
and Vasicek models. Tables [5] and [5] show the computed val ues of the callable bond in comparison to other 
metho d s. The columns BW, DF VL, BBKL, and F refe r to iBiittler and Waldvogel (|l996l) . Id'Halluin et~aT 



(|200ll) . iBen-Ameur et all (|2007l) . and Ide Frutosl (|2008h . Table [7] shows the results for the subordinated 



models in the present paper. We stress that neither of the alternative approaches in the literature is capable 
of handling jump- diffusion and pure jump models with state dependent jumps. The remarkable advantage of 
the eigenfunction expansion method is that it is entirely straightforward to move from pure diffusion models to 
jump- diffusion and pure jump models obtained by subordination by simply replacing the diffusion eigenvalues 
X n with subordinate eigenvalues <p(X n ). 

While Tables 1-7 provide results for the bond with the call option only to facilitate comparisons with the 
literature, Tables 8-10 provide the corresponding results for the bond that is both callable and putable. To 
generate this example of a bond with both options, we added the put option to the callable bond considered 
previously in this section. Table |S] gives the put prices we have assumed. Tabled presents results for break- 
even short rates for call and put options under the range of short rate models considered in this paper. Table 
[TU] presents the corresponding prices of the bond with both call and put options. 

6 Conclusion 

This paper proposed an efficient method to evaluate bonds with embedded options under a wide class of 
interest rate models, including the popular short rate diffusion models, as well as their time changed versions 
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with jumps. The method is based on the eigenfunction expansion of the pricing operator. Given the set 
of call and put dates, the callable and putable bond pricing function is the value function of a stochastic 
game with stopping times. Under some technical conditions, it is shown to have an eigenfunction expansion 
in eigenfunctions of the pricing operator with the expansion coefficients determined through a backward 
recursion. For CIR and Vasicek the method is orders of magnitude faster than the alternative approaches 
in the literature. In contrast to the alternative approaches in the literature that have so far been limited to 
diffusions, the method is equally applicable to short rate jump-diffusion and pure jump models constructed 
from diffusion models by Bochner's subordination with a Levy subordinator. In future work we plan to 
apply the eigenfunction expansion method of this paper to convertible bonds, where the stock price process 
is the stochastic variable driving the conversion and call decisions. 



7 Appendix 

Proof of Theorem 13.11 (i) For any / £ L 2 (I,m), Vtf € L 2 (I,m) for any t > 0. Since m is a finite 
measure on/, 1 € L 2 (/,m), so P{t,x) e L 2 (I,m) for any t > 0. Then by ((12) to ((II). C k ~ 1 (x) and F fc_1 (a;) 
are in L 2 (I,m). By induction using (TT51) . (|14D . and (TT5) . it can be shown that the value function V°(x) and 
the value functions V l (x) and the continuation value function C % (x) are in L 2 (I, m) for alii = k*, k — 1. 

(ii) The expressions for the continuation value function and the value function at the time of the bond 
issue are given in (TT5) and ([T5|) . By part (i), the eigenfunction expansion can be obtained from (JTJ) or ((5), 
where < +1 = (V t+1 ,<p n ), i = k*,...,k-l, and cf = (V k \<p n ) for n = 0,1,2,.... 

(iii) Under the conditions given in section 2.1 or 2.2, the eigenfunction expansion for the density given 
in © or (HH) holds. Then for / <E L 2 (7,m) and any x e I, 



p n OO 

'Ptfix) = f(y)p t (x,y)m(y)dy = f(y)y2er Xnt tp n (x)(p n {y)m(y)dy 

OO „ OO 

= e~ Kt Pn(x) / f(y)<Pn(y)m(y) dy = ^ /«e" A "VnW, 

n=0 J 1 ri=0 

where /„ = (/, ip n )- The interchange in the third equality is justified by the Dominated Convergence Theorem 
with the dominant function X^J^Lo e ~ Xnt \ l Pn(x)f(y)(p n (y)m(y)\: 



oo „ oo 

]T / e- A "* \<p n (x)f(y)<p n (y)m(y)\ dy < 

n=0^ 7 n=0 



-A„t 



Vn{x)\ ||/|U 2 H^n|U 2 



n=0 



The first inequality follows from the Cauchy- Schwartz inequality, and the last inequality follows from the 
bounds on eigenfunctions described in section 2.1 or 2.2 (e.g. |^ n (x)| < e A "*/ 2 y/pt(x, x) or |y> n (x)| < Ck) 
and the trace class condition ((2|) or ©. Hence, the eigenfunction expansion converges pointwise to Vtf(x) 
for any x G /. 

The eigenfunction expansions converge uniformly in x on compacts in / by the following: Let K be any 
compact subset of /. For x € K , 



-A„t 



n=M 



< E \fnC- Kt Vn(x)\ < \\fh* E e-^bnMI- 



n=M 



The last expression goes to as M goes to oo by the bounds on |y> n (a;)| and the trace class condition. 
Therefore, the eigenfunction expansions converge uniformly in x on compacts in I. 
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(iv) By (fT2|) . the values of c k are given by the coefficients of the eigenfunction expansion of the zero-coupon 
bond: 

c k n = (l + C) Pn , n = 0,1,2,... 

Suppose that we know the values of the coefficients c^ +1 , n = 0, 1, 2, .... Denote l( x ,y) as an indicator function 
that is 1 on the interval {x, y) and otherwise. 

4 = (V\^ n )+C Pn e~ x - s 

\m=0 / 

= K?p n {l, xf) +ljr c^e-^tfim, ¥>nl(xj,»f) J + *?Pn«> r) + Cp„ e - A " 5 

\m=0 / 

(OO OO \ 

J2 c^ 1 e" Am ' ! > m , ^ ^(if.x?)^ +Kfp n (x p i ,r) + C Pn e~ x " s 
m=0 m=0 / 

oo 

= K c iPn (l,xf) + Y, c^e-^w^ixlx?) +K?p n (x p i ,r) + C Pn e~ x " s . 



The last equality follows from (f,g) = J2™ =0 f n g n , where /„ = (f,tp n ) and g n = (g,ip n ) for f,g G 
I?{I,m). □ 

Proof of Proposition 13721 Suppose that (121)1) holds. From this we get that, with probability one, 

g-JXXiO^ds > e -/ t r(X 2 («))d* i 

Hence, the zero-coupon bond function P(t, x) is a positive decreasing function. 

We can also show that the continuation value function C l (x) and bond value function V 1 (x) are positive 
decreasing functions of x for all i = k*, k — 1. For i = k — 1, C k ~ 1 (x) is a constant multiple of a zero- 
coupon bond function by (IT21 and (fT5)) , so C 1 (a;) is a positive decreasing function. By (fH)) . the bond 
value function V^cc) is a positive decreasing function if C z (x) is a positive decreasing function since we 
already showed that zero-coupon bond functions are positive decreasing functions. It remains to show that 
for i = k* , k — 2 the continuation function C l (x) is a positive decreasing function given that V' l+1 (x) is a 
positive decreasing function. This is shown by the following. By (|26| , 



for JSfi(O) < X 2 (0). 

We show that there can be at most one solution x such that 

KP(S,x) = C l (x), 

where K is a constant (either Kf or Kf), for each i = k* , k — 1. It is first shown that there is at most 
one solution to the equation given by 

KP(8,x) = C k ~ x {x). 
Suppose that Xk-i is a solution to the equation. By (|13p. 



C k -\x) = (1 + C)E 
= (1 + C)E 



- f* r(X(s))ds 



-X*(rte— i) = x 



e-&K x ('»*'g(X(t k - 1 )) 



X(T k -l) = x 



(l + C)P(t fc _! - TTfc-x.a;)^**- 1 [ff(A:(t fc _i))|X(7)b_i) = x] , 
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where g{x) = E 



X(tk-i) = x is a decreasing function and E 1 * 1 - 1 denotes the expectation 

under the £fc_i-forward adjusted measure. Denote the probability measure under the -forward adjusted 
measure as Qk-i- Qk-i an d P are equivalent probability measures, so P[Xi(t) < X 2 {t) for all i > 0] = 1 
implies that <Q k _ x [X\{t) < X 2 (t) for all t > 0] = 1. Then E 1 "- 1 \g(X(t k -i))\X(Tk-i) = x] is a decreasing 
function of x. At x = x k —i, 

KP[t k -x - Tk-i,x) = (1 + C)P(t k -! - T k - 1 ,x)E tk -' [g(X(t k ^))\X{T k ^) = x] . 

Then K = (1 + C)E^ [g(X(t k ^))\X(T k -i) = z fe _i], so 

K > (1 + C)E tk -' [g(X(t k ^))\X(T k ^) = x] 

for all x > Xk-i- Then KP{t k -\ — r k -\,x) > C k ~ 1 (x) for all x > x k -\. 

It is shown next that there is at most one solution to the equation given by 

KP(S,x) = C l (x), 

for each i = k*, k — 2, where if is a constant (either Kf or Kf). Let Xi be such that KP(5, Xi) = C l (xi). 
By (T 



C l (a;) = £ 
= £ 



3-^* +V ^( s » ds ^ +1 (X(r i+1 ))|x(T i ) = 

r^ r W«w d ' fl (A:(t i ))x(r i ) = x 



= P{5,x)E t <[g{X{t i ))\X{ n )=x], 



where g(x) = E 



V l+1 (X (tj+i)) = cc is a decreasing function and E fi denotes the 



expectation under the ii-forward adjusted measure. Denote the probability measure under the ii-forward 
adjusted measure as Qj. Qi and P are equivalent probability measures, so P[Xl(£) < X2(t) for alH > 0] = 1 
implies that ®i[Xi(t) < X 2 (t) for all t > 0] = 1. Then E li [g(X{ti))\X{Ti) — x] is a decreasing function of 

X ■ -A_t X — X 'i ■ 

KP(6,x) = P{5 lX )E u [g(X(U))\X(n) = x] 

Then K = E u [g(X(ti))\X(n) = x l ],soK> E u [g(X(U))\X(n) = x] for all a; > Xi . Then KP(5,x) > C l {x) 
for all x > Xi. □ 



Proof of Proposition 14.11 For Laguerre polynomials, the forward shift property is 

— i (Q) - -L {a+1) (x) 
dx " ~ 1 j 



and the backward shift property is 

_d_ 
dx 

To derive the recursion for a^m(x), we first apply the backward shift, integrate by parts and then apply 
the forward shift. 



e- x x a L^(x)] = (n+^e-V-'L^ix). 



a n+l,m+l( X ) 



1 


TO + 


1 


1 




TO + 


1 


1 




TO + 


1 


1 





TO + 1 



e-yy^L^\y)d(L^ 1 (y) 
L { :l 1 {x)L^ +1) {x)e- x x a+1 + I e-yy^Lt +1 \x)L^\x)dy 



L^Ux)L^{x)e- x x^ + a^\x) 



(45) 
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For n > 1, m > 1, m ^ n, by solving the above equation for a"~^(x) and equating the expressions for 

< + ™( x ) and CW- we have 



m — n 

For n > 1, using the backward shift property, we have 



(lM(x)L%±1\x) - L${x)L^\x)) 



„(^) = ie-V+ 1 4t h 1 1) W- 



For n — m, we get the following from (1451) : 

a£> (a) = ~ [^(z^^K^ 1 + a^_i(^)] , (« > 1), a$ = 7(« + M> 

where j(a + 1, x) = J Q e~ v y a dy is the lower incomplete gamma function. 

The coefficients 6^ (x) are computed from equations (|46p and (|4T)) below. 
For n > 1, 



b&>(x) 



y a e- s yLl a \y)dy 
■ e -(-i)»I d [e-V +1 4- + i 1) (j/) 



s - 1 



(46) 



6^(x) 



^e-^4 a) (y)dy 
-Wa + Mx) 



(47) 
□ 



Proof of Proposition 14.21 For Hermite polynomials, the forward shift property is 

d 



d.r 



H n (x) = 2nH n ^i(x) 



and the backward shift property is 



_d_ 

d.r; 



e x H n (x) = -e x H n+1 (x). 



To derive the recursion for a n>m {x), we first apply the backward shift, integrate by parts and then apply 
the forward shift. 

dn+l,m+l(l) = 



H n+1 (y)H m+1 (y)e v dy 
H n+1 (y)d(e-y 2 H m (yY 

-OO 

= -H n+1 (y)H m {y)e-y 2 * +[ e~ y ' H m (y) d (H n+1 (y)) 
= -H n+1 (x)H m (x)e- x2 + f e- y2 H m (y)2{n + l)H n {y)dy 



-H n+1 (x)H m (x)e x +2(n + l)a rliTn (x) 



(48) 
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Exchanging the roles of n and m, we have 

2 

a m+ltn+ i(x) = -H m+1 {x)H n {x)e~ x + 2(m + \)a m . n {x). 

If rn 7^ n, then subtracting (|48|) from (|49|) and rearranging, we have 

, , _ g K (i)g m -nM - H m (x)H n+ i(x) _ x 2 
z(to — 71 ) 

Q-n,n(x) can be computed recursively as follows: 

a„, n {x) = -H n -i(x)H n (x)e~ x + 2na„_i !n _i(x), (n > 1), a , (x) = v^^V^x), 

where <j>(x) is the cumulative distribution function of a standard normal distribution. 
For n > 1, 



(49) 



6„(s,x) 



e su ' u H n {u)du 



e su d e~ u fln-i^ 

se s "-"~tf„_i(u)dw 

OO 

e* x_ar fl„_i(x) +sfe„_i(s,a;) 



b (s,x) = ieT^(Erf(i(2a ; - s)) + 1), 



where Erf (a;) is the error function. 



□ 
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Table 1: Call prices 



Exercise date Call Price 



hi = 


10.172 






1.025 






tl2 = 


11.172 






1.020 






*13 = 


12.172 






1.015 






tl4 = 


13.172 






1.010 






tl5 = 


14.172 






1.005 






<16 = 


10.172 to t 20 = 


19.172 




1.000 






Table 2: Parameter values 








Vasicek 




CIR 




K 






0.44178462 




0.14294371 




a 






0.13264223 




0.38757496 




e 






0.098397028 




0.133976855 




Table 3: Break-even short rates 


Time 


cm 


Vasicek 


SubCIR, JD 


SubCIR, PJ 


Sub Vasicek, JD 


Sub Vasicek, PJ 


T20 


0.03388791 


0.02706597 


0.03614163 


0.03672670 


0.03189678 


0.03348832 


T19 


0.01792789 


-0.01012520 


0.02292836 


0.02439808 


0.00299207 


0.00734621 


T18 


0.00978966 


-0.03655983 


0.01665424 


0.01758017 


-0.01809927 


-0.01208475 


T17 


0.00488209 


-0.05701483 


0.01161351 


0.01333251 


-0.03477951 


-0.02766935 


T16 


0.00157881 


-0.07350682 


0.00873978 


0.01047766 


-0.04847549 


-0.04061315 


T15 


n. a. 


-0.09100438 


n.a. 


n.a. 


-0.06370872 


-0.05539452 


T14 


n. a. 


-0.10481935 


n.a. 


n.a. 


-0.07568237 


-0.06698556 


T13 


n. a. 


-0.11653925 


n.a. 


n.a. 


-0.08590952 


-0.07694429 


T12 


n. a. 


-0.12671317 


n.a. 


n.a. 


-0.09485232 


-0.08570132 


Til 


n.a. 


-0.13566906 


n.a. 


n.a. 


-0.10277749 


-0.09350086 
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Table 4: Convergence and CPU times for pricing the callable bond with dynamic truncation level for initial 
short rate r = 0.05 



Pricing Error 


Average N at T20, m and to 


Maximum N at T20, rn and to 


CPU time (ms) 


cm 


1(T 


-5 


6.0, 3.4, 3.2, 3.1, 3.0, 4.0, 4.0, 4.0, 4.0, 4.0, 
2 


6, 5, 5, 5, 4, 4 


, 4, 4, 4, 4, 2 


1.1 


io- 


-6 


8.9, 7.0, 5.3, 5.4, 5.1, 5.0, 6.0, 6.0, 6.0, 6.0, 


9, 8, 8, 8, 7, 5 


, 6, 6, 6, 6, 2 


1.4 


10" 


-7 


2 

10.9, 11.0, 7.9, 9.0, 7.0, 7.0, 7.0, 7.0, 7.0, 3 


11, 11, 12, 11, 


11, 7, 7, 7, 7, 7, 3 


1.9 


Vasicek 


10" 


-5 


4.2, 5.8, 6.0, 4.3, 5.8, 5.9, 5.2, 5.2, 5.2, 5.2, 
2 


5, 6, 6, 6, 6, 6 


, 7, 7, 7, 7, 2 


0.8 


10- 


-6 


6.0, 8.3, 9.8, 9.8, 9.2, 9.0, 9.0, 9.3, 9.9, 10.0, 


6, 10, 10, 10, 


11, 10, 11, 11, 11, 11, 3 


1.3 


10" 


-7 


3 

6.1, 12.0, 12.0, 13.0, 13.0, 12.8, 12.9, 13.9, 


7, 13, 14, 14, 


14, 13, 13, 14, 14, 14, 3 


1.8 






13.8, 13.5, 3 








SubCIR, Jump-diffusion 


10" 


-S 


10.0, 9.9, 9.0, 5.9, 7.0, 7.0, 7.0, 6.0, 6.0, 6.0, 


10, 10, 11, 11, 


10, 7, 7, 6, 6, 6, 3 


2.1 


10" 


-6 


3 

11.9, 12.8, 13.4, 8.6, 9.7, 10., 8.0, 8.0, 8.0, 


12, 13, 14, 14, 


14, 10, 8, 8, 8, 8, 3 


2.7 






8.0, 3 








10" 


-7 


14.0, 15.8, 18.6, 19.2, 13.5, 16.0, 11.0, 10.0, 


14, 16, 19, 21, 


20, 16, 11, 10, 10, 10, 4 


4.0 






10.0, 10.0, 4 












SubCIR 


, Pure jump 






10" 


-5 


10.9, 10.8, 9.0, 5.7, 5.6, 5.0, 6.0, 6.0, 6.0, 


11, 11, 12, 12, 


11, 5, 6, 6, 6, 6, 3 


2.2 






6.0, 3 








10" 


-6 


12.9, 14.8, 18.3, 16.7, 10.9, 13.0, 9.0, 8.0, 


13, 15, 19, 20, 


20, 13, 9, 8, 8, 8, 3 


3.6 






8.0, 8.0, 3 








10" 


-7 


16.1, 27.7, 27.7, 35.0, 22.3, 31.0, 13.0, 12.0, 


17, 35, 28, 36, 


39, 31, 13, 12, 12, 12, 4 


8.7 






12.0, 12.0, 4 








SubVasicek, Jump-diffusion 


10" 


-5 


6.0, 6.3, 6.3, 7.9, 8.0, 7.3, 7.1, 7.2, 7.1, 7.1, 


6, 8, 10, 8, 8, 9, 8, 9, 9, 9, 3 


1.5 


10" 


-6 


3 

6.0, 12.0, 12.2, 13.2, 13.2, 13.2, 13.0, 13.9, 


7, 13, 14, 15, 15, 15, 15, 15, 15, 14, 3 


2.5 






14.0, 14.0, 3 








10" 


-7 


8.0, 18.3, 19.6, 20.6, 20.9, 21.0, 20.2, 20.0, 


8, 21, 22, 21, 23, 23, 22, 22, 22, 22, 4 


4.2 






21.6, 21.5, 4 








SubVasicek, Pure jump 


10" 


-5 


6.0, 10.3, 11.8, 12.3, 13.1, 13.0, 12.9, 12.6, 


6, 12, 14, 15, 16, 15, 15, 15, 15, 15, 3 


2.5 






13.7, 13.8, 3 








10" 


-6 


7.0, 22.3, 23.4, 25.5, 26.8, 27.8, 27.6, 27.6, 


8, 29, 30, 29, 29, 31, 32, 30, 30, 31, 4 


6.2 






28.7, 28.5, 4 








io- 


-7 


8.0, 44.6, 43.7, 47.8, 48.6, 47.2, 51.0, 51.1, 


8, 55, 52, 57, 53, 54, 55, 57, 57, 56, 5 


15.9 



54.0, 52.1, 5 
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Table 5: CIR model: Values of the callable bond for initial short rate r obtained by five methods 



r 


BW 


DFVL 




BBKL 


F 


this paper 


0.01 


0.9392 


0.93926 




0.93921 


0.93922 


0.939259 


0.02 


0.9159 


0.91598 




0.91595 


0.91596 


0.915992 


0.03 


0.8933 


0.89333 




0.89330 


0.89331 


0.893341 


0.04 


0.8712 


0.87127 




0.87125 


0.87125 


0.871290 


0.05 


0.8498 


0.84980 




0.84978 


0.84979 


0.849823 


0.06 


0.8289 


0.82890 




0.82888 


0.82889 


0.828923 


0.07 


0.8085 


0.80855 




0.80854 


0.80854 


0.808577 


0.08 


0.7887 


0.78874 




0.78873 


0.78873 


0.788769 


0.09 


0.7694 


0.76945 




0.76945 


0.76945 


0.769484 


0.10 


0.7507 


0.75067 




0.75067 


0.75067 


0.750708 


Table 6: Vasicek model: Values of the callable bond fo 


r initial short rate r obtained by four methods 


r 


BW 




DFVL 




BBKL 


this paper 


0.01 


0.8556 




0.84282 




0.84285 


0.842845 


0.02 


0.8338 




0.82627 




0.82630 


0.826294 


0.03 


0.8223 




0.81010 




0.81009 


0.810091 


0.04 


0.8062 




0.79420 




0.79423 


0.794230 


0.05 


0.7904 




0.77868 




0.77871 


0.778702 


0.06 


0.7749 




0.76348 




0.76351 


0.763502 


0.07 


0.7598 




0.74860 




0.74862 


0.748621 


0.08 


0.7450 




0.73403 




0.73406 


0.734053 


0.09 


0.7305 




0.71977 




0.71980 


0.719792 


0.10 


0.7163 




0.70578 




0.70583 


0.705830 


Table 7: 


Subordinated models: Values of the callable bond for 


initial short rate r 


obtained by the eigenfunc- 


tion expansion method 












r 


SubCIR, JD 




SubCIR, PJ 




Sub Vasicek, JD 


Sub Vasicek, PJ 


0.01 


0.967362 




0.972668 




0.874805 


0.884935 


0.02 


0.941069 




0.946130 




0.855193 


0.864408 


0.03 


0.915446 




0.920208 




0.835999 


0.844285 


0.04 


0.890481 




0.894892 




0.817216 


0.824562 


0.05 


0.866160 




0.870174 




0.798837 


0.805233 


0.06 


0.842470 




0.846044 




0.780854 


0.786293 


0.07 


0.819396 




0.822492 




0.763261 


0.767737 


0.08 


0.796927 




0.799510 




0.746050 


0.749559 


0.09 


0.775050 




0.777087 




0.729215 


0.731754 


0.10 


0.753752 




0.755215 




0.712749 


0.714318 



Table 8: Put prices 



Exercise date Put Price 



in = 10.172 


1.015 


tl2 = 11.172 


1.010 


tl3 = 12.172 


1.005 


tl4 = 13.172 


1.000 


tl5 = 14.172 


0.995 


t 16 = 10.172 to t 20 = 19.172 


0.990 
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Tabic 9: Break-even short rates for callable and putable bond 



Time 


cm 


Vasicck 


SubCIR, JD 


SubCIR, PJ 


SubVasicek, JD 


SubVasicek, PJ 


Call Option 


T20 


0.03388791 


0.02706597 


0.03614163 


0.03672670 


0.03189678 


0.03348832 


Tig 


0.03050674 


0.01653941 


0.03271682 


0.03328046 


0.02560234 


0.02738706 


TlS 


0.03032523 


0.01570707 


0.03248071 


0.03298851 


0.02523355 


0.02696967 


Tit 


0.03031566 


0.01565754 


0.03246480 


0.03296440 


0.02521294 


0.02694260 


Tie 


0.03031515 


0.01565469 


0.03246373 


0.03296242 


0.02521179 


0.02694084 


T15 


0.02494569 


0.01423308 


0.02715933 


0.02765770 


0.01905492 


0.02088314 


T14 


0.02447879 


0.01412248 


0.02662948 


0.02705660 


0.01851858 


0.02030185 


T13 


0.02427643 


0.01407361 


0.02641769 


0.02682992 


0.01830026 


0.02007824 


Til 


0.02409131 


0.01402853 


0.02623127 


0.02663907 


0.01810142 


0.01987946 


Til 


0.02390885 


0.01398410 


0.02604835 


0.02645309 


0.01790549 


0.01968408 


Put Option 


T20 


0.04534067 


0.04044891 


0.04765628 


0.04838597 


0.04477592 


0.04625085 


Tig 


0.04136813 


0.01957849 


0.04346118 


0.04402728 


0.03798709 


0.03955459 


T18 


0.04117866 


0.01857462 


0.04320875 


0.04371211 


0.03762279 


0.03914175 


T17 


0.04116872 


0.01851743 


0.04319187 


0.04368645 


0.03760252 


0.03911518 


T16 


0.04116820 


0.01851414 


0.04319074 


0.04368434 


0.03760139 


0.03911347 


T15 


0.03572256 


0.01708147 


0.03780731 


0.03830289 


0.03139350 


0.03301665 


T14 


0.03519281 


0.01694566 


0.03720163 


0.03761288 


0.3080348 


0.03238388 


T13 


0.03493847 


0.01688151 


0.03693765 


0.03733291 


0.03052750 


0.03210465 


T12 


0.03470234 


0.01682184 


0.03670090 


0.03709169 


0.03027123 


0.03185029 


Til 


0.03446938 


0.01676298 


0.03646824 


0.03685602 


0.03001840 


0.03159982 



Table 10: Values of the callable and putable bond for initial short rate r obtained by the eigenfunction 
expansion method 



r 


cm 


Vasicek 


SubCIR, JD 


SubCIR, PJ 


SubVasicek, 
JD 


SubVasicek, 
PJ 


0.01 


1.030391 


0.995407 


1.054194 


1.058549 


1.022068 


1.030678 


0.02 


1.004673 


0.975223 


1.025454 


1.029652 


0.998893 


1.006540 


0.03 


0.979637 


0.955474 


0.997443 


1.001420 


0.976211 


0.982876 


0.04 


0.955265 


0.936150 


0.970147 


0.973843 


0.954015 


0.959680 


0.05 


0.931540 


0.917242 


0.943553 


0.946911 


0.932295 


0.936946 


0.06 


0.908443 


0.898741 


0.917644 


0.920614 


0.911044 


0.914668 


0.07 


0.885958 


0.880639 


0.892409 


0.894942 


0.890253 


0.892840 


0.08 


0.864068 


0.862926 


0.867831 


0.869886 


0.869914 


0.871456 


0.09 


0.842758 


0.845594 


0.843898 


0.845435 


0.850019 


0.850510 


0.10 


0.822011 


0.828635 


0.820595 


0.821579 


0.830559 


0.829996 
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